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A Short Survey on Arithmetic Transforms 
and the Arithmetic Hartley Transform 
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Abstract 

Arithmetic complexity has a main role in the performance of algorithms for spectrum evaluation. Arithmetic transform 
theory offers a method for computing trigonometrical transforms with minimal number of multiplications. In this paper, 
the proposed algorithms for the arithmetic Fourier transform are surveyed. A new arithmetic transform for computing 
the discrete Hartley transform is introduced: the Arithmetic Hartley transform. The interpolation process is shown to 
be the key element of the arithmetic transform theory. 
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1 Introduction and Historical Background 

Despite the existence of fast algorithms for discrete trans¬ 
forms (e.g., fast Fourier transform, FFT), it is well known 
that the number of multiplications can significantly in¬ 
crease their computational (arithmetic) complexity. Even 
today, the multiplication operation consumes much more 
time than addition or subtraction. Table [I] brings the clock 
count of some mathematical operations as implemented for 
the Intel Pentium''’'^ processor. Observe that multiplica¬ 
tions and divisions can be by far more time demanding than 
additions, for instance. Sine and cosine function costs are 
also shown. 

This fact stimulated the research on discrete transform al¬ 
gorithms that minimize the number of multiplications. The 
Bhatnagar’s algorithm |laj . which uses Ramanujan num¬ 
bers to eliminate multiplications (however, the choice of 
the transform blocklength is rather limited), is an example. 
Parallel to this, approximation approaches, which perform a 
trade-off between accuracy and computational complexity, 
have been proposed [2aH4a] . 

Arithmetic transforms emerged in this framework as an 
algorithm for spectrum evaluation, aiming the elimination 
of multiplications. Thus, it may possess a lower compu¬ 
tational complexity. The theory of arithmetic transform 
is essentially based on Mobius function theorems (5^, of¬ 
fering only trivial multiplications, i.e., multiplications by 
{ — 1,0,1}. Therefore, only addition operations (except for 
multiplications by scale factors) are left to computation. 
Beyond the computational attractiveness, arithmetic trans¬ 
forms turned out to be naturally suited for parallel process¬ 
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Table 1: Clock count for some arithmetic instructions car- 
ried on a Pentium^^ processor [I] 


Operation 

Clock count 

add 

1-3 

sub 

1-3 

f add 

1-7 

f sub 

1-7 

mul (unsigned) 

10-11 

mul (signed) 

10-11 

div (unsigned) 

17-41 

div (signed) 

22-46 

fdiv 

39 

sin, cos 

17-137 


ing and VLSI design mis]. 

The very beginning ofthe research on arithmetic trans¬ 
forms dates back to 1903 when the German mathematician 
Ernst Heinrich BrunR depicted in Figure [I] published the 
Grundlinien des wissenschaftlichnen Rechnens [1], the sem¬ 
inal work in this field. In spite of that, the technique re¬ 
mained unnoticed for a long time. Forty-two years later, 
in Baltimore, U.S.A., the Hungarian Aurel Freidrich Wint- 
neJl, privately published a monograph entitled An Arith¬ 
metical Approach to Ordinary Fourier Series. This mono¬ 
graph presented an arithmetic method using Mobius func¬ 
tion to calculate the Fourier series of even periodic func¬ 
tions. 

After Wintner’s monograph, the theory entered again in 


^Bruns (1848-1919) earned his doctorate in 1871 under the super¬ 
vision of Weierstrass and Kummer. 

^A curious fact: Wintner was born in April 8th, 1903, in Budapest, 
the same year when Bruns had published his Grundlinien. Wintner 
died on January 15th, 1958, in Baltimore. 
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Figure 1: Ernst Heinrich Bruns: pioneer of arithmetic 
transform theory. Image in public domain. 

“hibernation” state. Not before 1988, Dr. Donald W. Tufts 
and Dr. Angaraih G. Sadasiv, independently, had rein¬ 
vented Wintner’s arithmetical procedure, reawaking the 
arithmetic transform. 

In the quest to implement it, two other researchers played 
an important role: Dr. Oved Shisha of the U.R.I. Depart¬ 
ment of Mathematics and Dr. Charles Rader of Lincoln 
Laboratories. They were aware of Wintner’s monograph 
and helped Tufts in many discussions. In 1988 The Arith¬ 
metic Fourier Transform by Tufts and Sadasiv was pub¬ 
lished in IEEE Acoustic, Speech, and Signal Processing 
(ASSP) Magazine [5]. 

Another breakthrough came in early 1990s is due 
to Emeritus Professor Dr. Irving S. Reed—the co¬ 
inventor of the widely used Reed-Muller (1954) and Reed- 
Solomon (1964) codes. Author of hundreds of publications. 
Dr. Reed made important contributions to the area of sig¬ 
nal processing. Specifically on arithmetic transforms, in 
1990 Reed, Tufts and co-workers provided two fundamental 
contributions [SIE]. In 0, a reformulated version of Tufts- 
Sadasiv approach, the arithmetic Fourier transform (AFT) 
algorithm was able to encompass a larger class of signals 
and to compute Fourier series of odd periodic functions as 
well as even periodic ones. The publication of the 1992 A 
VLSI Architecture for Simplified Arithmetic Fourier Trans¬ 
form Algorithm by Dr. Reed and collaborators in the IEEE 
Transactions on ASSP [6] was another crucial slash on the 
subject. Indeed, that paper was previously presented at 
the International Conference on Application Specific Array 
Processors held in Princeton. However, the 1992 publica¬ 
tion reached a vastly larger public, since it was published in 
a major journal. The new method, an enhancement of the 
last proposed algorithm [5] , was re-designed to have a more 
balanced and computationally efficient performance. As a 
matter of fact, Reed et al. proved that the newly proposed 
algorithm was identical to Bruns’ original method. 

When the AFT was introduced, some concerns on the 
feasibility of the AFT were pointed out [7] . The main issue 
dealt with the number of samples required by the algorithm. 
However, later studies showed that the use of interpolation 
techniques on a sub-sampled set (e.g., zero- and first-order 


interpolation) could overcome these difficulties [5]. 

The conversion of the standard 1-D AFT into 2-D ver¬ 
sions was just a matter of time. Many variants were pro¬ 
posed following the same guidelines of the 1-D case BM- 
Further research was carried out seeking different imple¬ 
mentations of the AFT. An alternative method [15] pro¬ 
posed a “Mobius-function-free AFT”. Iterative [16] and 
adaptative approaches m were also examined. In spite 
of that, the most popular presentations of the AFT are still 
those found in 016] • 

Although the main and original motivation of the arith¬ 
metic algorithm was the computation of the Fourier Trans¬ 
form, further generalizations were performed and the arith¬ 
metic approach was utilized to calculate other transforms. 
Dr. Luc Knockaert of Department of Information Technol¬ 
ogy at Ghent University, Belgium, amplified the Bruns pro¬ 
cedure, defining a generalized Mdbius transform [TSlIT^ . 
Moreover, four versions of the cosine transform was shaped 
in the arithmetic transform formalism |20j . 

Further generalization came in early 2000s with the defi¬ 
nition of the Arithmetic Hartley Transform (AHT) [2Tll22j . 
These works constituted an effort to make arithmetical pro¬ 
cedure applicable for the computation of trigonometrical 
transforms, other than Fourier transform. In particular 
the AHT computes the discrete Hartley transfornU: the 
real, symmetric, Fourier-like discrete transform defined in 
1983 by Emeritus Professor Ronald Newbold Bracewell in 
The Discrete Hartley Transform, an article published in the 
Journal of Optical Society of America. 

In 1988 and then the technological state-of-art was dra¬ 
matically different from that Bruns and Wintner found. 
Gomputational facilities and digital signal processing in¬ 
tegrated circuits made possible AFT to leave theoretical 
constructs and reach practical implementations. Since its 
inception in engineering, the AFT was recognized as tool 
to be implemented with VLSI techniques. Tufts himself 
had observed that AFT could be naturally implemented 
in VLSI architectures [2]. Implementations were proposed 
in [3l|6j|T2l|23]-|33]. Initial applications of the AFT took 
place in several areas: pattern matching techniques [34] . 
measurement and instrumentation [351136] . auxiliary tool 
for computation of z-transform [37l[38], and imaging [39] . 

This paper is organized in two parts. In Section B the 
mathematical evolution of the Arithmetic Fourier Trans¬ 
form is outlined. In Section (3] a summary of the major re¬ 
sults on the Arithmetic Hartley Transform is shown. Inter¬ 
polation issues are addressed and many points of the AFT 
are clarified, particularly the zero-order approximation. 


®Ralph Vinton Lyon Hartley (1888-1970) introduced his real inte¬ 
gral transform in a 1942 paper published in the Proceedings of I.R.E. 
The Hartley transform relates a pair of signals f(t) <—> F{v) by 

1 r°° 

F(v) = _ / /(i)(cos(i/t)-I-sin(i/t))di, 

\^27r J —oo 
1 

f{t) = _ / F{u){cos{iAt) sm{iAt))du. 

27r y _ oo 
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2 The Arithmetic Fourier Transform This is the finite version of the Mobius inversion for¬ 

mula [ 53 . A proof can be found in [5]. 

In this section, the three major developments of the arith¬ 
metic Fourier transform technique are presented. With em- 2 i Tufts-Sadasiv Approach 
phasis on the theoretical groundwork, the AFT algorithms 

devised by Tufts, Sadasiv, Reed et alii are briefly outlined. Consider a real even periodic function expressed by its 
In this work, ki\k 2 denotes that fci is a divisor of ^ 2 ; [-J is Fourier series, as seen below: 
the floor function and [•] is the nearest integer function. ^ 

v{t) ='^Vk{t). (7) 

/c=l 


Lemma 1 Let k, k' and m be integers. 

k-l 


m—O 


cos 2Trm- = 


k if k\k', 

0 otherwise 


and 


k-l 

E 

m=0 


sin I 27rm— = 0. 


Proof: Consider the expression Ylm=o 
k\k', yields 

k-l k-l 

T. = E 1 = 




m—0 

Otherwise, we have: 

k-l 


m—0 


(1) The components Vk{t) represent the harmonics of v{t), given 
by: 

Vk{t) = Ok ■ cos(27rA:t), (8) 

(2) where Ok is the amplitude of the fcth harmonic. 

It was assumed, without loss of generality, that v{t) had 
unitary period and null mean (oq = 0). Furthermore, con- 
When sider the N first harmonics as the only significant ones, in 
such a way that Vkit) = 0, for k > N (bandlimited ap¬ 
proximation). Thus the summation of Equation [7] might be 
constrained to N terms. 

Definition 2 The nth average is defined by 




m—0 

Therefore, it follows that: 

k-l 


m—0 


1 _ pi2irfe' 

,, = 0 . 
1 - e^'2-^ 


k if k\k', 

0 otherwise. 


n \ n J 


m—0 


for n = 1,2,..., A. Sn { t ) is null for n > N . 


(9) 


□ 


After an application of Equations [7] and [5] into IHl it 
yielded: 


By taking the real and imaginary parts, we conclude the 
proof. □ 

Definition 1 (Mobius /x-function) For a positive inte¬ 
ger n, 

(l ifn=l, 

fi{n) = < (—1)’' if n = rii^iPu Pi distinct primes, (3) 
lO if p^\n for some prime p. 

A relevant lemma related to /i-function is stated below. 

Lemma 2 

I ifn = l, 


n—1 


m 

27Tkt — 27rk — 
n 


E^(^) = 

d\n 


0 ifn>l. 


(4) 


m—0 

^ n—1 cxD 

=- EE®fc 

m—0 k—1 
^ oo n — 1 

= — Ok yy ^cos(27rfct) cos 

^ k—1 m—0 

sin(27r/ct) sin (^2'!rk —^ ^ 

1 °° 

= — 7 Ofe cos(27rfct) • 


27Tk 


m 


k^l 

oo 


Theorem 1 (Mobius Inversion Formula for Finite Series) 

Let n be a positive integer and fn a non-null seguence for 


= E 

n\k m=l 


In if n\k, I 
1 0 otherwise J 

lit), n = l,...,N. (10) 


I < n < A and null for n > N. If 

9 

then 


LAT/nJ 

~ y ' fkm 

k=l 

(5) 

IN/nj 

yy p.im)gmn- 

m—1 

(6) 


Proceeding that way, the nth average can be written in 
terms of the harmonics of p(t), instead of its samples (Defi- 
nition[7]). Since we assumed Vnit) = 0,n > N, only the first 
[A/nJ terms of Equation [TUI might possibly be nonnull. 

As a consequence the task was to invert Equation (TUI 
Doing so, the harmonics can be expressed in terms of the 
averages, S'„(t), which are derived from the samples of the 
signal p(t). The inversion is accomplished by invoking the 
□ Mobius inversion formula. 
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Theorem 2 The harmonics of v{t) can he obtained by: 


Vk{t) = fj,{m)Smkit), Vk=l,...,N. (11) 

m—1 

Proof: Some manipulation is needed. Substituting Equa¬ 
tion [ini in Equation [TTl it yields 

OO OO OO 

h{'ni)Smk{t) = Vkmn{t)- (12) 

m—1 m—1 n—1 

Now it is the tricky part of the proof. 

OO OO OO OO 

y] h{m)Y,vkmn(t) = y] y] H{m)vkmn{t) 

m—1 n—1 m—1 n—1 

\ (13) 

H{m) . 

According to Lemma [U the inner summation can only be 
null if j7 ^ = 1- In other words, the term Vk{t) is the 
only survivor of the outer summation and the proof is 
completed. □ 

The following aspects of the Reed-Tufts algorithm can be 
highlighted [5]: 

• This initial version of the AFT had a strong constraint: 
it can only handle even signals; 

• All computations are performed using only additions 
(except for few multiplications due to scaling); 

• The algorithm architecture is suitable for parallel pro¬ 
cessing, since each average is computed independently 
from the others; 




t=i 


A delay (shift) of aT in v{t) furnishes to the following ex¬ 
pression: 


j(t + aT) = y^ a„ cos f 27rn(^ -|- a)j + 

" ( t \ 

y^ bn sin f 27rn(— -f a) j 

n=l ^ 


n—1 

N 


n—1 


= y] c„(a) cos ( 27rn— j -|- 

n^l ^ ^ 

y^d„(a)sin f 27rn—j , 


where — 1 < a < 1 and 

Cri(a) = On cos(27rna) -I- bn sin(27rna), 
dn{of) = —an sin(27rna) -I- bn cos(27rna). 


(16) 


(17) 

(18) 


In the sequel, the computation of the Fourier coeffi¬ 
cients On and bn based on c„(a) is outlined. Meanwhile, 
the formula for the nth average (Tufts-Sadasiv) is updated 
by the next definition. 

Definition 3 The nth average is given by 

n — 1 

Snia)^-y2v(-T + aT), (19) 

n \ n J 

m=0 


where — 1 < a < 1. 

Now the quantities Cn{a) can be related to the averages, 
according to the following theorem. 


• The arithmetic transform theory is based on Fourier 
series, instead of the discrete Fourier transform. 

2.2 Reed-Tufts Approach 

Presented by Reed et al. in 1990 [5], this algorithm is a gen¬ 
eralization of Tuft-Sadasiv method. The main constraint of 
the latter procedure (handling only with even signals) was 
removed, opening path for the computation of all Fourier 
series coefficient of periodic functions. 

Let v{t) be a real T-periodic function, whose A^-term fi¬ 
nite Fourier series is given by 



Theorem 3 The coefficients Cn{a) are computed via 
Mobius inversion formula for finite series and are expressed 
by 

lN/n\ 

Cn(a) = ^ fJ,{l)Sln{a). (20) 

Proof: Substituting the result of Equation [!§] into Equa¬ 
tion M furnishes the following expression: 


N 


n — 1 


Sn{a) =y^Cfc(a)- y^ 


cos 


N 


m—0 
n — 1 


n 


sin 




m—O 


( 2'Kkm\ 

/ 2'Kkm\ 

\ n )' 


( 21 ) 


where Oq is the mean value of v{t). The even and odd 
coefficients of the Fourier series are a„ and respectively. 

Let v{t) denote the signal v{t) removed of its mean 
value ag. Consequently, 

v(t) = v(t) - ao 



A direct application of Lemma [T] yields 

[N/n] 

Sn(a) = y^ Cin(a). (22) 

1=1 


(15) 
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Invoking the Mobius inversion formula for finite series, the 
theorem is proved. □ 

Finally, the main result can be derived. 








Theorem 4 (Reed-Tufts) The Fourier series coefficients 
0,1 and bn are computed by 


— Cn(0), 

bn = {-iTCn 


f 1 


I 2^+2 


n = l,...,N, 


(23) 

(24) 


where k and m are determined by the factorization n = 

2'=(2m + 1). 

Proof: For a = 0, using Equation 1171 it is straightforward 
to show that a„ = c„(0). For a = ^ = 2^(2to-|-1), 

there are two sub-cases: m even or odd. 


For m = 2q, n = 2'^(4q + 1). Therefore, 


2'Kna = 2 tt 


2'=(4g-fl) ^ TT 

2^+2 - 27^9 + 2 • 


Consequently, substituting this quantity in Equa¬ 
tion [T71 yields 


1 


2k+2 


= a„ cos 

= bn- 


(2TTq 


+ 2 I + sin 


(2^9+1) 

(26) 


For m = 2q + 1, n = 2^(4g -f 3). It follows that 

2'=(4q-f3) 37r , , 

2'Kna = 27: -- = 27rg-f-^. (27) 

Again invoking the Equation 1171 the following expres¬ 
sion is derived. 


=anCOs(^2Kq+^^ + 

, , / 37r\ 

bn sin l2Kq+— j 

— bn- 


Joining these two sub-cases, it is easy to verify that 
bn = (-l)^Cn ' ^ 


2k+2 1 ■ 


Mr(N) = -N, 


2.3 Reed-Shih (Simplified AFT) 

Introduced by Reed et al [3] , this algorithm is an evolution 
of that one developed by Reed and Tufts. Surprisingly, in 
this new method, the averages are re-defined in accordance 
to the theory created by H. Bruns [J] in 1903. 

Definition 4 (Bruns Alternating Average) The 2nth 
Bruns alternating average, B 2 n(of), is defined by 

. 2n— 1 


m=0 


\2n 


(25) 


Invoking the definition of c„, applying Theorem [3] and 
Definition [3l the following theorem was derived. 

Theorem 5 The coefficients Cn(a) are given by the Mobius 
inversion formula for finite series as 


LwJ 

Cn(a) = ^ n(l) ■ B2nl(oi)- 

Z=l,3,... 


(33) 


(28) 


Proof: See [B]. □ 

Since a relation between the signal samples and the Bruns 
alternating averages was obtained, as well as an expression 
connecting the Bruns alternating averages to the c„ coeffi¬ 
cients, was available, few points are missing to compute the 
Fourier series coefficients. Actually, an expression relating 
the Fourier series coefficients (on and bn) to the coefficients 
Cn is sought. Examining Equation [T3 two conditions are 
distinguishable: 

• On = c„(0); 

• bn = Cn (^)- 

Calling Theorem [SI the next result was obtained. 

Theorem 6 (Reed-Shih) The Fourier series coefficients 
On and bn are computed by 


ao = 


T 


(29) 

□ 

The number of real multiplications and additions of this 
algorithm is given by [5] 


r 

v(t)dt, 

(34) 

LfJ 

E 

M(0^2ni(0), 

(35) 

1,3,5, 

L^J 

E 

1 .T R 

(A.)) , 

(36) 


(30) 


for n = 1 ,..., N. 

Proof: The proof is similar to the proof of Theorem H) □ 
For a blocklength N, the multiplicative and additive com¬ 
plexities are given by 


and 


Mr(N) = N, 


(37) 


Ar(N) = -Nf (31) 

respectively, where N is the blocklength of the transform. 


and 


Ar(N) = -Nf 


(38) 








respectively. 

The AFT algorithm proposed by Reed-Shih presents 
some improvements over previous algorithms: 

• The computation of both an and bn has roughly the 
same computational effort. The algorithm is more 
“balanced” than Reed-Tufts algorithm; 

• The algorithm is naturally suited to a parallel process¬ 
ing implementation; 

• It is computationally less complex that Reed-Tufts al¬ 
gorithm. 

2.4 An Example 

In this subsection, we draw some comments in connection to 
an example originally proposed by Reed et al. [3]. Let v{t) 
be a signal with period T = 1 s. Consider the computation 
of the Fourier series coefficients up to the 5th harmonic. 

According to the Reed-Shih algorithm, the coefficients an 
and bn of the Fourier series of v{t) are expressed by 


Oi 


' 10-10 - 1 ' 


■ ^ 2 ( 0 ) 

02 


0 10 0 0 


^4(0) 

03 

= 

0 0 10 0 


^ 6 ( 0 ) 

04 


0 0 0 1 0 


^ 8 ( 0 ) 

05_ 


0 0 0 0 1 


BioiO) 


and 


bi 


'10 10 - 1 ' 


■ 

62 


0 10 0 0 


Bi{\) 

bs 

= 

0 0 10 0 



64 


0 0 0 1 0 



65 . 


0 0 0 0 1 


_ Bw{^) _ 


Comparing these formulations with the ones of Reed-Tufts 
algorithm, one may note the balance in the computation 
of Qn and bn- Both coefficients are obtained through sim¬ 
ilar matrices. Table [2] relates the Bruns alternating aver¬ 
ages Bn(a) to their required samples. Notice that at least 
40 non-uniform time samples of v(t) are necessary to exactly 
compute the Bruns alternating averages. 

At this point, some observations are relevant: 

• This algorithm is not naturally suited for uniform sam¬ 
pling. 

• A uniform sampler utilized to obtain all the required 
samples would need a very high sampling rate. In the 
example illustrated here, a 120 Hz clock should be re¬ 
quired to sample the necessary points for the computa¬ 
tion of the Fourier series of a 1 Hz bandlimited signal. 

Certainly these observations appear to be disturbing and 
seems to jeopardize the feasibility of the whole procedure. 
However, it is important to stress that this procedure fur¬ 
nishes the exact computation of the Fourier series coeffi¬ 
cients. 

An empirical solution to circumvent this problem is to 
interpolate. An interpolation based on uniformly sampled 


Table 2: Necessary samples for the Bruns alternating aver¬ 
ages 


Bruns averages 

Sample time (s) 


^2(0) 

0,i 


^4(0) 



i?6(0) 

n 1 1 1 2 5 

’ 6’ 3’ 2’ 3’ 6 


BM 

r, 1 1 3 1 5 3 7 
8’ 4’ 8’ 2’ 8’ 4’ 8 


Sio(O) 

n 1 1 3 2 1 3 

10’ 5’ 10’ 5’ 2’ 5’ 

749 

10’ 5’ 10 

B2{\) 

1 3 

4’ 4 


Ba{\) 

13 5 7 

8’ 8’ 8’ 8 


Beih) 

1 1 5 7 3 11 

12 ’ 4 ’ 12 ’ 12 ’ 4 ’ 12 


B^ih) 

1 3 5 7 9 11 
16’ 16’ 16’ 16’ 16’ 16 

13 15 
’ 16’ 16 

Bw{^) 

1 3 1 7 9 11 

20’ 20’ 4’ 20’ 20’ 20’ 

13 3 17 19 
20’ 4’ 20’ 20 


points can be used to estimate the sample values required 
by AFT. Of course, this procedure inherently introduces 
computation errors. 

For example, assuming that the 1 Hz signal v{t) is sam¬ 
pled by a clock with period Tq = s. Hence, the following 
sample points are available: 



Table m shows, for example, that the computation of 
i? 4 ( 0 ) requires — among other samples — ^(i)) which 
is clearly not available. To overcome this difficulty, a 
rounding operation can be introduced. Thus, the sample 
V (^) can be used whenever the algorithm called v (|) 
([l0|;] /lO = 3/10). This rounding operation is also known 
as zero-order interpolation. 

The accuracy of the AFT algorithm is deeply associated 
with the sampling period Tq. If more precision is required, 
then one should expect to increase sampling rate, resulting 
in the introduction of smaller errors due to the interpola¬ 
tion scheme. Higher order of interpolation (e.g. first-order 
interpolation) can also be used to obtain more accurate es¬ 
timations of the Fourier series coefficients. The following 
trade-off is quite clear accuracy versus order of interpola¬ 
tion. 

However, for signals sampled at Nyquist rate (or close to), 
zero-order interpolation already leads to good results [5] . A 
detailed error analysis of interpolation schemes can be found 
in [5l [2011371140] . Further comments can be found in [3]. 

3 A New Arithmetic Transform 

Besides its numerical appropriateness [ 6 a], the discrete 
Hartley transform (DHT) has proved along the years to 
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be an important tool with several applications, such as 
biomedical image compression, OFDM/CDMA systems, 
and ADSL transceivers. Searching the literature, no men¬ 
tion about a possible “Arithmetic Hartley Transform” to 
compute the DHT was found. 

In this section, a condensation of the main results of 
the Arithmetic Hartley Transform is outlined. The method 
used to define the AHT turned out to furnish a new insight 
into the arithmetic transform. In particular, the role of in¬ 
terpolation is clarified. Additionally, it is mathematically 
shown that interpolation is a pivotal issue in arithmetic 
transforms. Indeed it determines the transform. 

A new approach to arithmetic transform is adopted. In¬ 
stead of considering uniformly sampled points extracted 
from a continuous signal v{t), the AHT is based on the 
purely discrete signal. Thus, the starting point of the de¬ 
velopment is the discrete transform definition, not the series 
expansion, as it was done in the development of the AFT al¬ 
gorithm. This approach is philosophically appealing, since 
in a final analysis a discrete transform relates two set of 
points, not continuous functions. 

Let V be an A^-dimensional vector with real ele¬ 
ments. The DHT establishes a pair denoted by v = 
[uo,-!;!,... o V = [Vb, lA,..., where the 

elements of the transformed vector V (i.e.. Hartley spec¬ 
trum) are defined by [5^ 



fc = 0,l,...,fV-I, (41) 


where casx = cosx -I- sin a; is Hartley’s “cosine and sine” 
kernel. The inverse discrete Hartley transform is then ( 6 ^ 

Uj = ^ 14-cas (-—) , i = 0, l,...,Af-l. (42) 

k=0 ^ ' 


An application of the inverse Hartley transform (Equa¬ 
tion |42) in Equation |44] offered: 


Sk 




cas 


2TTk'^ 

N 


(45) 


m—Ok'—O \ / 

Rearranging the summation order, simplifying, and calling 
Lemma [3l it yielded: 


Sk 


N-l k-1 

k'—O m—0 


cas I 27rm— 


L(iv-i)/fcj 

Ufc- 


s=0 


(46) 


For simplicity and without loss of generality, consider a 
signal V with zero mean value, i.e., ^ Vi = 0. Clearly, 

this consideration has no influence on the values of 14 , k ^ 
0. An application of the modified Mdbius inversion formula 
for finite series [5] is sufficient to obtain the final theorem 
to derive the Arithmetic Hartley Transform. According to 
Theorem [1] the result below follows. 

Theorem 7 (Reed et alii) If 

L(Af-l)/fcJ 

Sk= l<fc<A^-I, (47) 

S = 1 


then 

[(N-l)/k] 

Vk= 

1=1 

where ^{■) is Mdbius function. □ 

To illustrate the application of above theorem, consider 
an 8 -point DHT. Using Theorem[3 EauationHSl we obtain: 


Lemma 3 (Fundamental Property) The function 
cas(-) satisfies 


k-1 


m—0 


cas 2^m-- = 


k if k\k', 

0 otherwise. 


(43) 


Proof: It follows directly from Lemma [7J 


□ 


Similarly to the AFT theory, it is necessary to define av¬ 
erages Sk, calculated from the time-domain elements. The 
averages are computed by 

k — 1 

= fc = l,...,7V-l. (44) 

m—0 

This definition requires fractional index sampling. Analo¬ 
gously to the AFT methods, this fact seems to make further 
considerations impracticable, since only integer index sam¬ 
ples are available. This subtle question is to be treated in 
the sequel. For now, we assume that the fractional index 
sample are somehow available. 


Vi=Si-S2-S3-S5 + Se- Sr, 

V2 = S2-S4- Se, Us = As, 

V3 = A3 — Sq, Vq = Ag, 

Uj = A 4 , U 7 = A 7 . 

The component Ug = Ug can be computed directly from 
the given samples, since it represents the mean value of the 
signal Vb = I Figure [2] shows a diagram for this 

computation. 

The above theorem and equations completely specified 
how to compute the discrete Hartley spectrum. Addition¬ 
ally, the inverse transformation can also be established. The 
following result is straightforward. 

Corollary 1 Inverse discrete Hartley transform compo¬ 
nents can be computed by 

V{N-l)/^\ 

Vi= y] ti{l)crii, (49) 

1=1 

where a^ = = ^ 
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Figure 2: Diagram for computing the AHT for = 8. 
The boxes compute the averages and the multipliers imple¬ 
ment the scaling operation. The third layer accounts for 
the arithmetic computation based on Mbbius functions. 

The original Arithmetic Fourier Transform had identi¬ 
cal equations to those just derived for the Hartley trans¬ 
form (compare Equation 1101 and Equation 1471) . A question 
arises: since the equations are the same, which spectrum 
is actually being evaluated? Fourier or Hartley spectrum? 
A clear understanding of underlying arithmetic transform 
mechanisms will be possible in the next section. Once more 
the reader is asked to put this question aside for a while, 
allowing further developments to be derived. 

To sum it up, at this point two major questions are ac¬ 
cumulated: (i) How to handle with fractional indexes? and 
(ii) How does same formulae result in different spectra? In¬ 
terestingly, both questions had the same answer. 

The arithmetic transform algorithm can be summarized 
in four major steps: 

1. Index generation, i.e., calculating the indexes of neces¬ 
sary samples 

2. Fractional index samples handling, which requires in¬ 
terpolation; 

3. Computation of averages: Sk = ^ Sm=o 

4. Computation of spectrum by Mobius Inversion For¬ 
mula: 14 = ^{l)Ski. 


In the rest of this paper, the step two is addressed. In the 
sequel, a mathematical method, explaining the importance 
of the interpolation process in the arithmetic algorithms, is 
derived. 

3.1 Interpolation 

Arithmetic transform theory usually prescribes zero- or 
first-order interpolation for the computation of spectral ap¬ 
proximations 01113H]- In this section, it is shown that 
an interpolation process based on the known components 
(integer index samples) characterizes the definition of the 
fractional index components, Vr, r 0 N. This analysis al¬ 
lows a more encompassing perception of the interpolation 
mechanisms and gives mathematical tools for establishing 
validation constraints to such interpolation process. In ad¬ 
dition, brief comments on the trade-off between accuracy 
and computational cost required by interpolation process 
close the section. 

3.1.1 Ideal Interpolation 

What does a fractional index discrete signal component re¬ 
ally mean? The value of Vr for a noninteger value r, r ^ N, 
can be computed by 



Defining the Hartley weighting function by 



the value of the signal at fractional indexes can be found 
utilizing an A^-order interpolation expressed by: 

N-l 

i^O 

It is clear that each transform kernel can be associated to 
a different weighting function. Consequently, a different in¬ 
terpolation process for each weighting function is required. 
In the arithmetic transform formalism, the difference from 
one transform to another resides in its interpolation process. 

The weighting functions satisfy 'Wi{r) = 1. If r 

is an integer number, then the orthogonality properties of 
cas(-) function make Wr{r) = 1 and Wi{r) = 0 (Vi r). 
Therefore, no interpolation is needed. 

After some trigonometrical manipulation, the interpola¬ 
tion weights for several kernels is expressed by closed for¬ 
mulae. Let Sa(-) be the sampling function, 

Sa(a:) = 


sin(ai)/a;, x 0, 
1, a; = 0. 


(53) 
























































r=0,0.1 ...15 



Figure 3: Hartley weighting functions used to interpolate 
uio.i and U 10.5 {N = 32 blocklength). 


Proposition 1 The N-point discrete Fourier cosine, 
Fourier sine, and Hartley transforms have interpolation 
weighting functions given by: 


Fourier cosine Kernel 


Wi 


,-(r) - J- + J i 

*V' ; — 27V ^ N 12 


Sa( '^ W^^ 27r(z-r)) 


2 Sa{7r{i—r)/N) 


Fourier sine kernel 

N- 1/2 f 1 Sa(K^f^27r(i-r)) 
N ] 2 Sa(7r(i-r)/N) 


W, 


'i(r) = 


Hartley kernel 

1 I N- 1/2 


= 27v + 


N 


Sa(7r('i—r) /N) 


2N 


cot 




1 Sa( ^ ;y^^^ 27r(i+r)) ) 

2 Sa(7r(i+r)/A^) J 


1 Sa( ^ ;^^^^ 27r(i+r)) '[ 

2 Sa{n{i-\-r)/N) 


/■ 


1 cos(i^-j^27r(i+r)) 
2N sin{7T /N) 

□ 




(b) Hartley kernel 


To exemplify, Figure [3] shows two weighting functions 
used to compute uio.i and wio.s in the arithmetic Hartley 
transform. Figure 0] shows the weighting profile for the 
Fourier cosine and Hartley kernels. These functions are 
calculated by closed formulae. 

With this proposition, the mathematical description of 
the AHT algorithm is completed. The derived formulae 
furnish the exact value of the spectral components. On 
the other hand, the computational complexity of the ideal 
interpolation implementation is similar to the direct imple¬ 
mentation, i.e., computing the transform by its plain defi¬ 
nition: Vk = cas To address this issue, 

a non-ideal interpolation scheme is proposed. 

3.1.2 Non-Ideal Interpolation 

According to the index generation (m^), the number R 
of points that require interpolation is upper bounded by 
R < d — 1. This sum represents the number of sam¬ 

ples with fractional index. Consequently, this approach can 
be attractive for large non-prime blocklength N with great 
number of factors, because it would require a smaller num¬ 
ber of interpolations. 


Figure 4: These curve families represent the weighting 
profile Wi{r) for the Fourier cosine kernel (a) and Hartley 
kernel (b). The parameters A^ = 16, r = 0, 0.1... 15, and 
i = 0... 15 are employed. Observe that i G N and r G K. 
The maximum values are achieved ati = 0ori = A^/2 = 8 
(central peak) and they corresponds to the unity. The value 
of the local maxima are exactly 1 / 2 , occurring when r is 
integer. Note that each curve is almost null everywhere, 
except at the vicinities of i ~ r or i ss A^ — r. 

The next task is to find simpler formulae for the weighting 
functions, assuming large blocklength. Instead of using the 
exact weighting functions, the limit N^oo is examined and 
used to derive asymptotic approximations of the weighting 
function. 

Proposition 2 A continuous approximation for the inter¬ 
polation weighting function for sufficiently large N is given 
by: 

Fourier cosine kernel 

Wi{r) « _|_ Sa(2TT-q-|-7-)) 

Fourier sine kernel 
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Wt{r) 

Hartley kernel 

Wi{r) '■ 


^ Sa(27r(i—r)) Sa(27r(2+r)) 

' 2 ~ 2 

Sa(27r(*-r)) + ig^. 


□ 

In terms of the Hilbert transform and the Sa(-) func¬ 
tion, the asymptotic weighting function for Hartley kernel 
is given by 

Wi(r) ss Sa(27r(i — r)) — T-Lil | Sa(27r(I -|- r)) |, (54) 

or alternatively, 

Wi{r) K, Sa(27r(i — r)) — Ca(27r(z -|- r))— (55) 

T-Lil^5{2{r -\- i)) 

where T-Lil denotes the Hilbert transform, Ca(a;) = , for 

a; ^ 0, is the co-sampling function and S{x) is the Dirac 
impulse. 


Zero-order Interpolation. Rounding the fractional 
index provides the zero-order interpolation. The estimated 
(interpolated) signal Vj is then expressed by Vj = T[j], 
where [•] is a function which rounds off its argument to 
its nearest integer. Examining the asymptotic behavior of 
the weighting function for Fourier cosine kernel, we derive 
the following results: 


w*(r)w0, Vi 7 ^ [r],IV - [r], 
Sa (27r([r] — r)) I 

2 2 ’ 

Sa (27r([r] — r)) I 

2 2 ' 


W[r] (r) 
WN-[r]{r) 


(56) 


Under the above assumptions, the Equation [52] furnishes: 


Vr W W[r](r)u[r] -k WW-[r](l’)l^Ar-[r'] 

1 1 (57) 

Thus, for even signal (vk = VN-k) , the approximated 
value of the interpolated sample is roughly given by Vr ~ 
vy\- It is straightforward to see that the influence of odd 
part of the signal vanishes in the zero-order interpolation. 
Besides zero-order interpolation is “blind” to the odd com¬ 
ponent of a signal. In fact, as show by the set of Equa¬ 
tions EH zero-order interpolation is an (indeed good) ap¬ 
proximation to the cosine asymptotic weighting function. 
This puts some light on the role of the zero-order inter¬ 
polations and its relation to the earliest versions of the 
arithmetic Fourier transform, which can not analyze odd 
periodic signals. 

Zero-order interpolation has been intuitively employed in 
previous work by Tufts, Reed et alii WM- Hsu, in his 
Ph.D. dissertation, derives an analysis of first-order inter¬ 
polation effect [38) . 



k 


Figure 5: The discrete Hartley transform computed by def¬ 
inition (solid line) and by arithmetic transform algorithm 
(dotted line). Simulation data: f{t) = cos(907rt) (t — 
t = 0...I, A^ = 32. 


Interpolation Order. Let Mm be a set with the m 
(< N) most significant coefficients Wi{r). For zero-order 
interpolation, m = I. Increasing the value of m, the in¬ 
terpolation process would gradually be improved, because 
more coefficients would be retained. Proceeding in this way, 
the following calculation performs a non-ideal interpolation: 


Vr 


- ^ w,{r) ■ Vi, 
ri 

‘ iGMm 


(58) 


where rj = ^ normalization factor. 

Figure El presents a 32-point discrete Hartley transform 
of a particular signal computed by the plain definition of 
the DHT and by the arithmetic transform method using 
m = 2. 


4 Conclusions 

This paper supplied a short survey on arithmetic trans¬ 
forms. The arithmetic Fourier transform is explained and 
its three main 1-D versions were described. Furthermore, 
some comments on implementation, challenging points, and 
advantages of the arithmetic transforms were discussed via 
simple examples. 

In this paper, the introduction of the AHT emphasized 
the key point of the arithmetic transforms: the interpola¬ 
tion process. It was shown that the fundamental equations 
of the arithmetic transform algorithms were essentially the 
same (regardless the kernel). This property could open path 
to the implementation of “universal transformers”. In this 
type of construction, the circuitry would be the same for 
several transforms, except the interpolation module. A dif¬ 
ferent interpolation module would reflect different trans¬ 
form (Fourier, Hartley, Fourier cosine). Finally, this paper 
could be taken as a starting point to those who want to 
investigate arithmetic transforms. 
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